#include "HY1C_out.h"
#include <sys/types.h>
#include <unistd.h>

/* ============================================================================ */
/* sst.c - functions for retrieval of sea surface temperature                   */
/*                                                                              */
/* Written By: B. Franz, NASA/SIMBIOS, August 2003.                             */
/*                                                                              */
/* ============================================================================ */

#include "l12_proto.h"
#include "hdf_utils.h"

static float sstbad = BAD_FLT;
static int32 evalmask = 0;

typedef struct stat_struct {
  float min;
  float max;
  float avg;
  int   cnt;
} statstr;

/* sensor-specific error table (SSES) */
#define NSSTMAX   8
#define NDAYMAX   2
#define NQUARMAX  4
#define NSENZMAX  4
#define NDIFFMAX  4
#define NLATMAX   5
#define NQUALMAX  4
#define NSSESDIM  7

typedef struct ssestab_struct {

    int nsst;
    int nday;
    int nquar;
    int nsenz;
    int ndiff;
    int nlat;
    int nqual;

    float sst [NSSTMAX ];
    float senz[NSENZMAX];
    float diff[NDIFFMAX];
    float lat [NLATMAX ];
    float qual[NQUALMAX];

    float bias[NQUALMAX][NLATMAX ][NDIFFMAX][NSENZMAX][NQUARMAX][NDAYMAX][NSSTMAX ];
    float stdv[NQUALMAX][NLATMAX ][NDIFFMAX][NSENZMAX][NQUARMAX][NDAYMAX][NSSTMAX ];

} ssestabstr;

static ssestabstr sses_sst;
static ssestabstr sses_sst4;

/* scans of computed quantities */
static float  *sst        = NULL;
static float  *sst4       = NULL;
static int16  *flags_sst  = NULL;
static int16  *flags_sst4 = NULL;
static int8   *qual_sst   = NULL;
static int8   *qual_sst4  = NULL;
static float  *bias_sst   = NULL;
static float  *bias_sst4  = NULL;
static float  *stdv_sst   = NULL;
static float  *stdv_sst4  = NULL;

/* precomputed indices */
static int32_t   recnumSST   = -1;
static int    initSST     =  0;
static int    haveSST4    =  0;
static int    haveSST     =  0;
static int    haveSSES    =  1;
static int    haveRed     =  0;
static int    ib39        = -1;
static int    ib40        = -1;
static int    ib11        = -1;
static int    ib12        = -1;
static int    ibred       = -1;
static int    nbvis       = -1;
static float  satred;

/* quality test parameters */
static int   cldbox     =    3;
static float cldthresh  = 0.01;
static int   btbox      =    3;
static float hisenz     = 55.0;
static float vhisenz    = 75.0;
static float Btmin      = -4.0;
static float Btmax      = 33.0;
static float SSTmin     = -2.0;
static float SSTmax     = 45.0;
static float dBtmin     =  0.0;
static float dBtmax     =  3.6;
static float dBt4min    =  0.0;
static float dBt4max    =  8.0;
static float SSTdiff    =  3.0;
static float SSTvdiff   =  5.0;
static float SST4diff1  =  0.8;
static float SST4diff2  =  1.0;
static float Bt11unif1  =  0.7;
static float Bt12unif1  =  0.7;
static float Bt11unif2  =  1.2;
static float Bt12unif2  =  1.2;
static float Bt39unif1  =  0.7;
static float Bt40unif1  =  0.7;
static float Bt39unif2  =  1.2;
static float Bt40unif2  =  1.2;
static float dBtrefmin  = -1.1;
static float dBtrefmax  = 10.0;

static float Bt11devavg = 0.10;   /* test not in use */
static float Bt12devavg = 0.15;
static float Bt39devavg = 0.10;
static float Bt40devavg = 0.15;

/* flag bit settings */
#define SSTF_ISMASKED    0x0001           /* Pixel aready masked          */
#define SSTF_BTBAD       0x0002           /* Brights are bad              */
#define SSTF_BTRANGE     0x0004           /* Brights are out-of-range     */
#define SSTF_BTDIFF      0x0008           /* Brights are too different    */
#define SSTF_SSTRANGE    0x0010           /* SST outside valid range      */
#define SSTF_SSTREFDIFF  0x0020           /* Different from reference     */
#define SSTF_SST4DIFF    0x0040           /* Different from SST4          */
#define SSTF_SST4VDIFF   0x0080           /* Very different from SST4     */
#define SSTF_BTNONUNIF   0x0100           /* BT window nonuniform         */
#define SSTF_BTVNONUNIF  0x0200           /* BT window very nonuniform    */
#define SSTF_BT4REFDIFF  0x0400           /* 4um BT diff relative to ref  */
#define SSTF_REDNONUNIF  0x0800           /* Red band spatial nonuniform  */
#define SSTF_HISENZ      0x1000           /* Sensor zenith high           */
#define SSTF_VHISENZ     0x2000           /* Sensor zenith very high      */
#define SSTF_SSTREFVDIFF 0x4000           /* Very diff from reference     */
#define SSTF_SPARE2      0x8000           /* Spare                        */

static float latwin = 2.5; /* half the overlap for lat band coeffs */

static int StartOfMonth[2][12] = {{0,31,59,90,120,151,181,212,243,273,304,334},
				  {0,31,60,91,121,152,182,213,244,274,305,335}};

static int pdbg1 = 360;
static int ldbg1 = 1760;
static int pdbg2 = 800;
static int ldbg2 = 245;

/* ---------------------------------------------------------------------------------------- */
/* load_sses_sst() - loads the specified SSES (sensor-specific error stats) table           */
/* ---------------------------------------------------------------------------------------- */
void load_sses_sst(int32_t sensorID, char *file, ssestabstr *sses)
{
    int32 sd_id;
    int32 sds_id; 
    int32 status;
    int32 rank; 
    int32 nt; 
    int32 nattrs;
    int32 dims [NSSESDIM]; 
    int32 start[NSSESDIM]; 

    char  name   [H4_MAX_NC_NAME]  = "";
    char  sdsname[H4_MAX_NC_NAME]  = "";

    memset(start,0,NSSESDIM*sizeof(int32));

    if (strcmp(file,"") == 0) {
        HY1C_out("\nNo SSES data provided for this sensor.\n");
        haveSSES = 0;
        return;
    }

    HY1C_out("\nLoading SSES table from %s\n",file);

    /* open the file and initiate the SD interface */
    sd_id = SDstart(file, DFACC_RDONLY);
    if (sd_id == -1) {
        HY1C_out("-E- %s:  Error opening file %s.\n", __FILE__,file);
        exit(1);
    }

    /* read the bias and standard deviation */

    strcpy(sdsname,"bias");
    sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
    status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
    if (rank != NSSESDIM) {
        HY1C_out("-E- %s: Table dimensions for %s do not match expectation.\n",__FILE__,sdsname);
        exit(1);
    }
    status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->bias);
    if (status != 0) {
        HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
        exit(1);
    }
    status = SDendaccess(sds_id);

    strcpy(sdsname,"stdv");
    sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
    status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
    if (rank != NSSESDIM) {
        HY1C_out("-E- %s: Table dimensions for %s do not match expectation.\n",__FILE__,sdsname);
        exit(1);
    }
    status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->stdv);
    if (status != 0) {
        HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
        exit(1);
    }
    status = SDendaccess(sds_id);


    /* save the table dimensions */
    sses->nqual = dims[0];
    sses->nlat  = dims[1];
    sses->ndiff = dims[2];
    sses->nsenz = dims[3];
    sses->nquar = dims[4];
    sses->nday  = dims[5];
    sses->nsst  = dims[6];


    /* read the table indice ranges */

    strcpy(sdsname,"sst");
    sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
    status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
    status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->sst);
    if (status != 0) {
        HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
        exit(1);
    }
    status = SDendaccess(sds_id);
    sses->nsst = dims[0];

    strcpy(sdsname,"senz");
    sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
    status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
    status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->senz);
    if (status != 0) {
        HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
        exit(1);
    }
    status = SDendaccess(sds_id);
    sses->nsenz = dims[0];

    strcpy(sdsname,"BTdiff");
    sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
    status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
    status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->diff);
    if (status != 0) {
        HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
        exit(1);
    }
    status = SDendaccess(sds_id);
    sses->ndiff = dims[0];

    strcpy(sdsname,"lat");
    sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
    status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
    status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) sses->lat);
    if (status != 0) {
        HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
        exit(1);
    }
    status = SDendaccess(sds_id);
    sses->nlat = dims[0];

    /* terminate access to the SD interface and close the file */
    status = SDend(sd_id);
}


/* ----------------------------------------------------------------------------------- */
/* comp_sses_sst() - computes sensor-specific error stats for SST                      */
/* ----------------------------------------------------------------------------------- */
void comp_sses_sst(ssestabstr *sses, float diff, float sst, float solz, float senz, 
               int doy, float lat, int iqual, float *bias, float *stdv)
{
    int i, iday, isst, idiff, isenz, iquar, ilat;

    /* get table indices */
    iquar = doy/91;
    if (solz > 90.0)
        iday = 1;
    else
        iday = 0;
    if ((evalmask & SSTMODS) != 0) iday = 1;   // use night data for all cases (7/22/2009)
    for (i=1; i<sses->nlat; i++) 
        if (lat < sses->lat[i])
	  break;
    ilat = i-1;
    for (i=1; i<sses->nsenz; i++) 
        if (senz < sses->senz[i])
	  break;
    isenz = i-1;
    for (i=1; i<sses->nsst; i++) 
        if (sst < sses->sst[i])
	  break;
    isst = i-1;
    for (i=1; i<sses->ndiff; i++) 
        if (diff < sses->diff[i])
	  break;
    idiff = i-1;

    /* set bias and standard deviation */
    if (iqual < sses->nqual) {
        *bias = sses->bias[iqual][ilat][idiff][isenz][iquar][iday][isst];
        *stdv = sses->stdv[iqual][ilat][idiff][isenz][iquar][iday][isst];
        if ((evalmask & SSTMODS) != 0) *bias += 0.17;  // add skin-temp offset (7/22/2009)
    } else {
        *bias = -999.0;
        *stdv = -999.0;
    }

    /*
    HY1C_out("%f %f %f %f\n",senz,lat,sst,diff);
    HY1C_out("%d %d %d %d %d %d %d\n",iqual,ilat,idiff,isenz,iquar,iday,isst);
    HY1C_out("%f %f\n",*bias,*stdv);
    */
}


/* ----------------------------------------------------------------------------------- */
/* btrefdiff() - scan-dependent test on 4um BT temps to reference (RSMAS)              */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
float btrefdiff(int32_t ip, float BT22, float BT23)
{
    static float dpix   = 13.0;
    static int   ntab   = 105;
    static float tref[] = {
                4.44, 4.22, 4.02, 3.97, 3.90, 3.73, 3.55, 3.44, 3.35, 3.29,
		3.18, 3.09, 3.03, 2.99, 2.93, 2.88, 2.81, 2.77, 2.74, 2.70,
		2.67, 2.64, 2.58, 2.56, 2.54, 2.51, 2.46, 2.46, 2.43, 2.41,
		2.38, 2.35, 2.35, 2.33, 2.32, 2.29, 2.27, 2.25, 2.24, 2.26,
		2.24, 2.21, 2.20, 2.21, 2.19, 2.19, 2.18, 2.17, 2.17, 2.16,
		2.17, 2.15, 2.16, 2.16, 2.15, 2.16, 2.16, 2.16, 2.16, 2.16,
		2.17, 2.15, 2.14, 2.14, 2.17, 2.21, 2.20, 2.19, 2.19, 2.19,
		2.18, 2.21, 2.26, 2.30, 2.31, 2.30, 2.34, 2.32, 2.35, 2.40,
		2.42, 2.45, 2.42, 2.50, 2.55, 2.58, 2.65, 2.71, 2.77, 2.82,
		2.89, 2.96, 3.02, 3.10, 3.16, 3.22, 3.31, 3.38, 3.45, 3.55,
		3.66, 3.80, 3.94, 4.08, 4.30};

    int32_t  itab = (int) ip/dpix;
    float wt   = ip/dpix - itab;
    float diff;

    if (itab < ntab-1)
      diff = BT22 - BT23 - (tref[itab]*(1.0-wt) + tref[itab+1]*wt);
    else
      diff = BT22 - BT23 - tref[itab];

    /* abs() to match miami code, but I think it's wrong
    return(fabs(diff));
    */
    return(diff);  /* fixed for collect 5 */
}

/* ----------------------------------------------------------------------------------- */
/* set_flags_sst4() - set quality flags for long wave sea surface temperature           */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
void set_sses_sst (l2str *l2rec)
{
    static int firstCall = 1;

    int32_t  npix = l2rec->npix;
    float sstref;
    int32_t  ip, ipb;
    float Bt11;
    float Bt12;
    float senz;
    float solz;
    float lat;
    float dBt;
    int   doy;

    /* set each pixels bias and stdev (assumes sst and qual computed) */
    for (ip=0; ip<npix; ip++) {
        ipb  = ip*NBANDSIR;
        senz = l2rec->senz[ip];
        solz = l2rec->solz[ip];
        lat  = l2rec->lat [ip];
        doy  = *l2rec->day;
        Bt11 = l2rec->Bt[ipb+ib11];
        Bt12 = l2rec->Bt[ipb+ib12];
        dBt  = Bt11-Bt12;
        comp_sses_sst(&sses_sst,dBt,sst[ip],solz,senz,doy,lat,qual_sst[ip], 
		      &bias_sst[ip], &stdv_sst[ip]);
    }
}


/* ----------------------------------------------------------------------------------- */
/* set_flags_sst4() - set quality flags for long wave sea surface temperature           */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
void set_sses_sst4(l2str *l2rec)
{
    static int firstCall = 1;

    int32_t  npix = l2rec->npix;
    float sstref;
    int32_t  ip, ipb;
    float Bt39;
    float Bt40;
    float senz;
    float solz;
    float lat;
    float dBt;
    int   doy;

    /* set each pixels bias and stdev (assumes sst and qual computed) */
    for (ip=0; ip<npix; ip++) {
        ipb  = ip*NBANDSIR;
        senz = l2rec->senz[ip];
        solz = l2rec->solz[ip];
        lat  = l2rec->lat [ip];
        doy  = *l2rec->day;
        Bt39 = l2rec->Bt[ipb+ib39];
        Bt40 = l2rec->Bt[ipb+ib40];
        dBt  = btrefdiff(ip,Bt39,Bt40);
        comp_sses_sst(&sses_sst4,dBt,sst4[ip],solz,senz,doy,lat,qual_sst4[ip], 
		      &bias_sst4[ip], &stdv_sst4[ip]);
    }
}


/* ----------------------------------------------------------------------------------- */
/* init_sst() - Initialize for SST processing                                          */
/* B. Franz, SAIC, August 2005                                                         */
/* ----------------------------------------------------------------------------------- */
void init_sst(l2str *l2rec)
{
    int32_t npix = l2rec->npix;
    int  i;

    evalmask = l2rec->input->evalmask;

    nbvis = NBANDS;
    ibred = bindex_get(  678);
    ib39  = bindex_get( 3959) - l2rec->nbands;
    ib40  = bindex_get( 4050) - l2rec->nbands;
    ib11  = bindex_get(11000) - l2rec->nbands;
    ib12  = bindex_get(12000) - l2rec->nbands;

    if (ib11 < 0 || ib12 < 0) 
        haveSST = 0;
    else
        haveSST = 1; 

    if (ib39 < 0 || ib40 < 0) 
        haveSST4 = 0;
    else
        haveSST4 = 1; 

    if (!haveSST && !haveSST4) {
        fHY1C_out(stderr,"-E- %s line %d: no SST bands found.\n",
        __FILE__,__LINE__);
        exit(1);
    }

    /* this is the value that a saturated radiance would be assigned */
    if (ibred < 0) 
        haveRed = 0;
    else {
        haveRed = 1;
	//        satred = 1000.0/l2rec->input->gain[ibred]-1.0;
        satred = 1000.0 - 1.0;
    }

    /* allocate private arrays for a single scan line */
    if (haveSST) {
        sst  = (float*) calloc(npix,sizeof(float));
        flags_sst  = (int16*) calloc(npix,sizeof(int16));
        qual_sst   = (int8*) calloc(npix,sizeof(int8));
        bias_sst   = (float*) calloc(npix,sizeof(float));
        stdv_sst   = (float*) calloc(npix,sizeof(float));
        load_sses_sst(l2rec->sensorID,l2rec->input->sstssesfile,&sses_sst );
    }

    if (haveSST4) {
        sst4 = (float*) calloc(npix,sizeof(float));
        flags_sst4 = (int16*) calloc(npix,sizeof(int16));
        qual_sst4  = (int8*) calloc(npix,sizeof(int8));
        bias_sst4  = (float*) calloc(npix,sizeof(float));
        stdv_sst4  = (float*) calloc(npix,sizeof(float));
        load_sses_sst(l2rec->sensorID,l2rec->input->sst4ssesfile,&sses_sst4);
    }

    initSST = 1;
}



/* ----------------------------------------------------------------------------------- */
/* sst_ran() - Determine if we have already run sst for this scan line                 */
/* ----------------------------------------------------------------------------------- */
int sst_ran(int recnum)
{
    if ( recnum == recnumSST )
        return 1;
    else
        return 0;
}


/* ----------------------------------------------------------------------------------- */
/* read_sst_mirror() - reads sst correction for MODIS/Terra mirror-side cooling (RSMAS)*/
/* B. Franz, SAIC, 11 August 2005                                                      */
/* ----------------------------------------------------------------------------------- */
void read_sst_mirror(int32 sensorID, char *filename, int year, int day, float *dsst)
{
    char mission[4];

    FILE *fp;
    char line [80];
    char odate [8];
    char sdate [8];
    int  found = 0;
    float lastdsst;

    *dsst = 0.0;

    switch (sensorID) {
      case MODIST: 
      case HMODIST: 
        strcpy(mission,"TER");
        break;
      default:
        return;
    }

    /* Open the file */
    if ( (fp = fopen(filename,"r")) == NULL ) {
        fHY1C_out(stderr,
                 "-E- %s line %d: unable to open %s for reading\n",
                __FILE__,__LINE__,filename);
        exit(1);
    }

    /* Form date string */
    sprintf(odate,"%4d%03d",year,day);

    /* Loop through to find closest time */
    while ( fgets( line, 80, fp ) ) {
        if ( strncmp(line,mission,3) == 0 ) {
	    lastdsst = *dsst;
            sscanf(line,"%4s %7s %f",mission,sdate,dsst);
            found = 1; /* we'll use the last available */
            if (strcmp(odate,sdate) <= 0)
                break;
	}
    }

    fclose(fp);

    HY1C_out("\nLoading mirror-side SST coefficients from %s.\n",filename);
    HY1C_out("Using mirror-side SST correction of %f from %s.\n",*dsst,sdate);

    if (found == 0) {
        fHY1C_out(stderr,
                 "-E- %s line %d: unable to locate mirror-side SST correction for %s in %s\n",
                __FILE__,__LINE__,odate,filename);
        exit(1);
    }

    return;
}


/* ----------------------------------------------------------------------------------- */
/* read_sst_coeff() - reads sst coefficients for the specified date                    */
/*                                                                                     */
/* B. Franz, SAIC, 11 August 2005                                                      */
/* ----------------------------------------------------------------------------------- */
void read_sst_coeff(int32 sensorID, char *filename, int32_t year, int32_t day, 
     float bounds[6][2], float coef[6][4])
{
    char mission[5] = "";

    FILE *fp;
    char line [80] = "";
    char odate [8] = "";
    char sdate [8] = "";
    char edate [8] = "";
    int32  found = 0;
    int32  indx  = 0;
    int32  tmonth = 0;
    int32  month;
    int32  leap;

    switch (sensorID) {
      case MODISA: 
      case HMODISA: 
        strcpy(mission,"AQUA");
        break;
      case MODIST: 
      case HMODIST: 
        strcpy(mission,"TERR");
        break;
      default:
        strcpy(mission,"AQUA");
        break;
    }

    /* Open the file */
    if ( (fp = fopen(filename,"r")) == NULL ) {
        fHY1C_out(stderr,
                 "-E- %s line %d: unable to open %s for reading\n",
                __FILE__,__LINE__,filename);
        exit(1);
    }


    if ((evalmask & SSTMODS) != 0) {

	/* Find month */
	leap = (isleap(year) == TRUE? 1: 0);
	for (tmonth=11; tmonth >= 0; tmonth--) {
	    if (day >= StartOfMonth[leap][tmonth]) {
		break;
	    }
	}

	/* Loop through to find 6 sets of coeffs for this month */
	indx=0;
	while ( fgets( line, 80, fp ) ) {
	    if ( strncmp(line,mission,4) == 0 ) {
		sscanf(line,"%4s %i %f %f %f %f %f %f",
		mission,&month,&bounds[indx][0],&bounds[indx][1],
		&coef[indx][0],&coef[indx][1],&coef[indx][2],&coef[indx][3]);
		if (month == tmonth+1) {
		    indx++;
		    if (indx == 6) {
			found = 1;
			break;
		    }
		}
	    }
	}

    } else {

	/* Form date string */

	uselastyear:
	sprintf(odate,"%4d%03d",year,day);


	/* Loop through to find bounding times, for 2 sets of coeffs */

	indx=0;
	while ( fgets( line, 80, fp ) ) {
	    if ( strncmp(line,mission,4) == 0 ) {
		sscanf(line,"%4s %7s %7s %f %f %f %f",
		mission,sdate,edate,&coef[indx][0],&coef[indx][1],&coef[indx][2],&coef[indx][3]);
		if ((strcmp(odate,sdate) >= 0 && strcmp(odate,edate) <= 0) || strcmp(edate,"0000000") == 0) { 
		    indx++;
		    if (indx == 2) {
			found = 1;
			break;
		    }
		}
	    }
	}

	if (found == 0 && year > 2004) {
	    HY1C_out("Warning: No SST coefficients available for %s, reverting to previous year.\n",odate);
	    year--;
	    rewind(fp);
	    goto uselastyear;
	}
    }


    fclose(fp);
    found=1;
    if (found) {

	if ((evalmask & SSTMODS) != 0) {

	    HY1C_out("Loading SST lat band coefficients from %s:\n",filename);
	    for (indx=0; indx<6; indx++) {
		HY1C_out("%i %6.1f %6.1f %6.3f %6.3f %6.3f %6.3f\n",month,bounds[indx][0],bounds[indx][1],
			coef[indx][0],coef[indx][1],coef[indx][2],coef[indx][3]);
	    }

	} else {
	    HY1C_out("Loading SST coefficients from %s:\n",filename);
	    HY1C_out("%s %s %6.3f %6.3f %6.3f %6.3f\n",sdate,edate,coef[0][0],coef[0][1],coef[0][2],coef[0][3]);
	    HY1C_out("%s %s %6.3f %6.3f %6.3f %6.3f\n\n",sdate,edate,coef[1][0],coef[1][1],coef[1][2],coef[1][3]);
	}

    } else {
        fHY1C_out(stderr,
                 "-E- %s line %d: unable to locate valid SST coefficients for %s in %s\n",
                __FILE__,__LINE__,odate,filename);
        exit(1);
    }

    return;
}


/* ----------------------------------------------------------------------------------- */
/* read_sst4_coeff() - reads 4um sst coefficients for the specified date               */
/*                                                                                     */
/* B. Franz, SAIC, 11 August 2005                                                      */
/* ----------------------------------------------------------------------------------- */
void read_sst4_coeff(int32 sensorID, char *filename, int32_t year, int32_t day, 
     float bounds[6][2], float coef[6][5])
{
    char mission[5];

    FILE *fp;
    char line [80];
    char odate [8];
    char sdate [8];
    char edate [8];
    int  found = 0;
    int  tmonth;
    int  month;
    int  indx;
    int32  leap;

    switch (sensorID) {
      case MODISA: 
      case HMODISA: 
        strcpy(mission,"AQUA");
        break;
      case MODIST: 
      case HMODIST: 
        strcpy(mission,"TERR");
        break;
      default:
        strcpy(mission,"AQUA");
        break;
    }

    /* Open the file */
    if ( (fp = fopen(filename,"r")) == NULL ) {
        fHY1C_out(stderr,
                 "-E- %s line %d: unable to open %s for reading\n",
                __FILE__,__LINE__,filename);
        exit(1);
    }


    if ((evalmask & SSTMODS) != 0) {

	sprintf(odate,"%4d%03d",year,day);

	/* Find month */
	leap = (isleap(year) == TRUE? 1: 0);
	for (tmonth=11; tmonth >= 0; tmonth--) {
	    if (day >= StartOfMonth[leap][tmonth]) {
		break;
	    }
	}

	/* Loop through to find 6 sets of coeffs for this month */
	indx=0;
	fHY1C_out(stderr, "  looking for month %d mission %s\n",tmonth+1,mission);
	while ( fgets( line, 80, fp ) ) {
	    if ( strncmp(line,mission,4) == 0 ) {
		sscanf(line,"%4s %i %f %f %f %f %f %f %f",
		mission,&month,&bounds[indx][0],&bounds[indx][1],
		&coef[indx][0],&coef[indx][1],&coef[indx][2],&coef[indx][3],&coef[indx][4]);
		if (month == tmonth+1) {
		    indx++;
		    if (indx == 6) {
			found = 1;
			break;
		    }
		}
	    }
	}

	fclose(fp);

//    } else if (sensorID == MODIST || sensorID == HMODIST) {

    } else {

	/* Form date string */

	uselastyear:
	sprintf(odate,"%4d%03d",year,day);

	/* Loop through to find bounding times */

	while ( fgets( line, 80, fp ) ) {
	    if ( strncmp(line,mission,4) == 0 ) {
		sscanf(line,"%4s %7s %7s %f %f %f %f",mission,sdate,edate,&coef[0][0],&coef[0][1],&coef[0][2],&coef[0][3]);
		if ((strcmp(odate,sdate) >= 0 && strcmp(odate,edate) <= 0) || strcmp(edate,"0000000") == 0) {
		    found = 1;
		    break;
            }
	    }
	}

	if (found == 0 && year > 2004) {
	    HY1C_out("Warning: No SST4 coefficients available for %s, reverting to previous year.\n",odate);
	    year--;
	    rewind(fp);
	    goto uselastyear;
	}

	fclose(fp);

//    } else {
//
//	/* use these sst4 coeffs instead of coeff file for AQUA */
//
//	coef[0][0] = 0.81;
//	coef[0][1] = 1.03;
//	coef[0][2] = 0.39;
//	coef[0][3] = 0.90;
//	coef[0][4] = 0.21;
//
//	found = 1;
    }
    found = 1;
    if (found) {

	if ((evalmask & SSTMODS) != 0) {
	    HY1C_out("Loading SST4 lat band coefficients from %s:\n",filename);
	    for (indx=0; indx<6; indx++) {
		HY1C_out("%i %6.1f %6.1f %6.3f %6.3f %6.3f %6.3f %6.3f\n",month,bounds[indx][0],bounds[indx][1],
			coef[indx][0],coef[indx][1],coef[indx][2],coef[indx][3],coef[indx][4]);
	    }
	} else {
	    HY1C_out("Loading SST4 coefficients from %s:\n",filename);
	    HY1C_out("%s %s %6.3f %6.3f %6.3f %6.3f %6.3f\n",sdate,edate,coef[0][0],coef[0][1],coef[0][2],coef[0][3],coef[0][4]);
	}
    } else {
        fHY1C_out(stderr,
                 "-E- %s line %d: unable to locate valid SST4 coefficients for %s in %s\n",
                __FILE__,__LINE__,odate,filename);
        exit(1);
    }

    return;
}


/* ----------------------------------------------------------------------------------- */
/* sstcloud() - uses red band to test homogeneity in nx x ny box (RSMAS)               */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
int sstcloud(int32_t ip, int nx, int ny, float thresh)
{
    extern l1qstr l1que;
    extern l1str *my_l1_data;
    int32_t   nscan  = 99999;
//    int32_t   nscan  = l1que.nq;
    int32_t   npix   = l1que.r[0].npix;
    int    is     = nscan/2;
    int32_t   ip1, ip2;
    int32_t   is1, is2;
    int32_t   i, j;
    int32_t   ipb;
    float  rhot;
    int    flag = 0;
    int    cnt  = 0;
    float  maxv = -9999.;
    float  minv =  9999.;

    if (!haveRed)
        return(flag);

    /* make sure code is not inconsistent NQMIN in l12_parms.h */
    if (nscan < ny) {
        HY1C_out("-E- %s line %d: L1 queue size of %d is too small for requested homogeneity test of %d x %d.\n", 
               __FILE__,__LINE__,l1que.nq,nx,ny);
        exit(1);
    }

    /* algorithm is only valid for daytime */
    if (l1que.r[is].solz[ip] >= SOLZNIGHT)
        return(flag);

    /* compute pix/scan limits for the ROI */
    is1 = MIN(MAX(0,is-ny/2),nscan-1);
    is2 = MAX(MIN(nscan-1,is+ny/2),0);
    ip1 = MIN(MAX(0,ip-nx/2),npix-1);
    ip2 = MAX(MIN(npix-1,ip+nx/2),0);


    /* compute max and min normalized reflectance in window */
    for (j=is1; j<=is2; j++) for (i=ip1; i<=ip2; i++) {
        ipb = i*nbvis + ibred;
	//HY1C_out("%d %d %f %f\n",i,j,l1que.r[j].Lt[ipb],satred);
        if (l1que.r[j].Lt[ipb] > 0.0 && l1que.r[j].Lt[ipb] < satred) {
            rhot = PI * l1que.r[j].Lt[ipb] / l1que.r[j].Fo[ibred]
                 / l1que.r[j].tg_sol[ipb]
                 / l1que.r[j].tg_sen[ipb]
                 / l1que.r[j].t_sol   [ipb]
                 / l1que.r[j].t_sen   [ipb]
	         / cos(l1que.r[j].solz[ip]/RADEG);
	    maxv = MAX(maxv,rhot);
	    minv = MIN(minv,rhot);
            cnt++;
        }
    }

    /* added for collect 5 (if all saturated but 1, then cloud) */
    if (cnt < 2 || (maxv-minv) > thresh) {
        flag = 1;
    }

   return(flag);
}


/* ----------------------------------------------------------------------------------- */
/* btboxstats() - uses red band to test homogeneity in nx x ny box (RSMAS)             */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
int32_t btboxstats(int32_t ip, int ib, int nbands, int32_t nx, int32_t ny, statstr *statrec)
{
    extern l1qstr l1que;

    int32_t   nscan  = l1que.nq;
    int32_t   npix   = l1que.r[0].npix;
    int32_t   ip1, ip2;
    int32_t   is1, is2;
    int32_t   is, i, j;
    int32_t   ipb;
    float  Bt;


    /* make sure code is not inconsistent NQMIN in l12_parms.h */
    if (nscan < ny) {
        HY1C_out("-E- %s line %d: L1 queue size of %d is too small for  %d x %d box stats.\n", 
               __FILE__,__LINE__,l1que.nq,nx,ny);
        exit(1);
    }

    /* Compute queue scan limits for the ROI */
    is  = nscan/2;
    is1 = MIN(MAX(0,is-ny/2),nscan-1);
    is2 = MAX(MIN(nscan-1,is+ny/2),0);

    /* compute pixel neighbor limits for the ROI */
    ip1 = MIN(MAX(0,ip-nx/2),npix-1);
    ip2 = MAX(MIN(npix-1,ip+nx/2),0);

    /* initialize output rec */
    statrec->min =  9999.;
    statrec->max = -9999.;       
    statrec->avg = 0.0;
    statrec->cnt = 0;

    /* compute stats for window */
    for (j=is1; j<=is2; j++) for (i=ip1; i<=ip2; i++) {
        ipb = i*nbands + ib;
        Bt  = l1que.r[j].Bt[ipb];
        if (Bt > 0.0) {
	    statrec->max = MAX(statrec->max,Bt);
	    statrec->min = MIN(statrec->min,Bt);
	    statrec->avg += Bt;
	    statrec->cnt++;
        }
    }
    if (statrec->cnt > 0)
        statrec->avg /= statrec->cnt;
        
    return(statrec->cnt);
}



/* ----------------------------------------------------------------------------------- */
/* sstmasked() - returns 1 if pixel was already masked (SST processing skipped)        */
/* ----------------------------------------------------------------------------------- */
int sstmasked(l2str *l2rec, int32_t ip)
{
    if ( (l2rec->flags[ip] & LAND)    != 0 ||
         (l2rec->flags[ip] & NAVFAIL) != 0 )
/*
         l2rec->pixnum[ip] < 8             ||
         l2rec->pixnum[ip] > 1345)
*/
        return(1);
    else 
        return(0);
}


/* ----------------------------------------------------------------------------------- */
/* set_flags_sst() - set quality flags for long wave sea surface temperature           */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
void set_flags_sst(l2str *l2rec)
{
    static int firstCall = 1;

    int32_t  npix = l2rec->npix;
    float sstref;
    int32_t  ip, ipb;
    float Bt11;
    float Bt12;
    float senz;
    float solz;
    float dBt;
    float dSST;

    statstr statrec;

    /* check each pixel in scan */
    for (ip=0; ip<npix; ip++) {

        /* SST not processed */
        if (sstmasked(l2rec,ip)) {
	    flags_sst[ip] |= SSTF_ISMASKED;
            continue;
	}

        ipb = ip*NBANDSIR;
        senz = l2rec->senz[ip];
        solz = l2rec->solz[ip];
        Bt11 = l2rec->Bt[ipb+ib11];
        Bt12 = l2rec->Bt[ipb+ib12];

        /* BT cound not be computed (radiance out-of-range) */
        if (Bt11 < BT_LO+0.1 || Bt11 > BT_HI-0.1 || 
            Bt12 < BT_LO+0.1 || Bt12 > BT_HI-0.1) {
            flags_sst[ip] |= SSTF_BTBAD;
            continue;
        }

        /* check BT range */
        if (Bt11 < Btmin || Bt11 > Btmax || 
            Bt12 < Btmin || Bt12 > Btmax)
            flags_sst[ip] |= SSTF_BTRANGE;

        /* check BT diff */
        dBt = Bt11-Bt12;
        if (dBt < dBtmin || dBt > dBtmax)
            flags_sst[ip] |= SSTF_BTDIFF;

        /* check SST range */
        if (sst[ip] < SSTmin || sst[ip] > SSTmax)
            flags_sst[ip] |= SSTF_SSTRANGE;

        /* check SST difference with references */
        dSST = sst[ip]-l2rec->sstref[ip];
	if ((evalmask & SSTMODS) != 0) {
          /* evaluate change to cold-test only */
          if (dSST < -SSTdiff)
            flags_sst[ip] |= SSTF_SSTREFDIFF;
          if (dSST < -SSTvdiff)
            flags_sst[ip] |= SSTF_SSTREFVDIFF;
	} else {
          if (l2rec->solz[ip] >= SOLZNIGHT) {
            if (fabs(dSST) > SSTdiff)
                flags_sst[ip] |= SSTF_SSTREFDIFF;
	  } else {
            if (dSST < -SSTdiff || dSST > (SSTdiff+1) )
                flags_sst[ip] |= SSTF_SSTREFDIFF;
	  }
          if (fabs(dSST) > SSTvdiff)
            flags_sst[ip] |= SSTF_SSTREFVDIFF;
	}

        /* check SST difference with 4um SST */
        /* set BT4REFDIFF flag based on sst4 flags */
        if (haveSST4 && sst4[ip] != sstbad) {
            dSST = fabs(sst[ip]-sst4[ip]);
            if (dSST > SST4diff1) 
                flags_sst[ip] |= SSTF_SST4DIFF;
            if (dSST > SST4diff2) 
                flags_sst[ip] |= SSTF_SST4VDIFF;
            flags_sst[ip] |= (flags_sst4[ip] & SSTF_BT4REFDIFF);
	}

        /* check sensor zenith limits */
        if (senz > hisenz)
            flags_sst[ip] |= SSTF_HISENZ;
        if (senz > vhisenz)
            flags_sst[ip] |= SSTF_VHISENZ;

        /* if SST is cold, check for clouds (collect 5) */
        if (sst[ip]-l2rec->sstref[ip] <= -1.0)
	    if (sstcloud(ip,cldbox,cldbox,cldthresh) == 1)
	        flags_sst[ip] |= SSTF_REDNONUNIF;

	/* check homogeneity of BT */
        if (btboxstats(ip,ib11,NBANDSIR,btbox,btbox,&statrec) > 0) {
	    if ((statrec.max - statrec.min) > Bt11unif1)
	        flags_sst[ip] |= SSTF_BTNONUNIF;
	    if ((statrec.max - statrec.min) > Bt11unif2)
	        flags_sst[ip] |= SSTF_BTVNONUNIF;
	    /*
            if (fabs(Bt11 - statrec.avg) > Bt11devavg)
	        flags_sst[ip] |= SSTF_BTDEVAVG;
	    */
	}
        if (btboxstats(ip,ib12,NBANDSIR,btbox,btbox,&statrec) > 0) {
	    if ( (statrec.max - statrec.min) > Bt12unif1)
	        flags_sst[ip] |= SSTF_BTNONUNIF;
	    if ( (statrec.max - statrec.min) > Bt12unif2)
	        flags_sst[ip] |= SSTF_BTVNONUNIF;
	    /*
            if (fabs(Bt12 - statrec.avg) > Bt12devavg)
	        flags_sst[ip] |= SSTF_BTDEVAVG;
	    */
	}
                
    }
    return;
}

/* ----------------------------------------------------------------------------------- */
/* set_flags_sst4() - set quality flags for long wave sea surface temperature           */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
void set_flags_sst4(l2str *l2rec)
{
    static int firstCall = 1;

    int32_t  npix = l2rec->npix;
    float sstref;
    int32_t  ip, ipb;
    float Bt39;
    float Bt40;
    float senz;
    float solz;
    float dBt;
    float dSST;

    statstr statrec;

    /* check each pixel in scan */
    for (ip=0; ip<npix; ip++) {

        /* SST not processed */
        if (sstmasked(l2rec,ip)) {
	    flags_sst4[ip] |= SSTF_ISMASKED;
            continue;
	}

        // ipb = ip*NBANDSIR;
        // 20180822
        ipb = ip*2;
        senz = l2rec->senz[ip];
        solz = l2rec->solz[ip];
        Bt39 = l2rec->Bt[ipb+ib39];
        Bt40 = l2rec->Bt[ipb+ib40];

        /* if BT cound not be computed (radiance out-of-range) */
        if (Bt39 < BT_LO+0.1 || Bt39 > BT_HI-0.1 || 
            Bt40 < BT_LO+0.1 || Bt40 > BT_HI-0.1) {
            flags_sst4[ip] |= SSTF_BTBAD;
            continue;
        }

        /* check BT range */
        if ( Bt39 < Btmin || Bt39 > Btmax || 
             Bt40 < Btmin || Bt40 > Btmax )
            flags_sst4[ip] |= SSTF_BTRANGE;

        /* check BT diff */
        dBt = Bt39-Bt40;
        if (dBt < dBt4min || dBt > dBt4max)
            flags_sst4[ip] |= SSTF_BTDIFF;

        /* check BT diff against BT reference */
        dBt = btrefdiff(ip,Bt39,Bt40);
        if (dBt < dBtrefmin || dBt > dBtrefmax)
            flags_sst4[ip] |= SSTF_BT4REFDIFF;

        /* check SST range */
        if (sst4[ip] < SSTmin || sst4[ip] > SSTmax)
            flags_sst4[ip] |= SSTF_SSTRANGE;

        /* check SST difference with reference */
        dSST = sst4[ip]-l2rec->sstref[ip];
	if ((evalmask & SSTMODS) != 0) {
          /* evaluate change to cold-test only */
          if (dSST < -SSTdiff)
            flags_sst4[ip] |= SSTF_SSTREFDIFF;
          if (dSST < -SSTvdiff)
            flags_sst4[ip] |= SSTF_SSTREFVDIFF;

	} else {
          if (l2rec->solz[ip] >= SOLZNIGHT) {
            if (fabs(dSST) > SSTdiff)
                flags_sst4[ip] |= SSTF_SSTREFDIFF;
	  } else {
            if (dSST < -SSTdiff || dSST > (SSTdiff+1) )
                flags_sst4[ip] |= SSTF_SSTREFDIFF;
	  }
          if (fabs(dSST) > SSTvdiff)
            flags_sst4[ip] |= SSTF_SSTREFVDIFF;
	}
    
        /* check SST4 difference with 11um SST */
        if (sst[ip] != sstbad) {
	    dSST = fabs(sst[ip]-sst4[ip]); 
            if (dSST > SST4diff1) 
                flags_sst4[ip] |= SSTF_SST4DIFF;
            if (dSST > SST4diff2) 
                flags_sst4[ip] |= SSTF_SST4VDIFF;
	}

        /* check sensor zenith limits */
        if (senz > hisenz)
            flags_sst4[ip] |= SSTF_HISENZ;
        if (senz > vhisenz)
            flags_sst4[ip] |= SSTF_VHISENZ;

        /* if 11um SST is cold, check for clouds (collect 5) */
        if (sst[ip] != sstbad && sst[ip]-l2rec->sstref[ip] <= -1.0)
	    if (sstcloud(ip,cldbox,cldbox,cldthresh) == 1)
	        flags_sst4[ip] |= SSTF_REDNONUNIF;

	/* check homogeneity of BT */
        if (btboxstats(ip,ib39,NBANDSIR,btbox,btbox,&statrec) > 0) {
	    if ((statrec.max - statrec.min) > Bt39unif1)
	        flags_sst4[ip] |= SSTF_BTNONUNIF;
	    if ((statrec.max - statrec.min) > Bt39unif2)
	        flags_sst4[ip] |= SSTF_BTVNONUNIF;
	    /*
            if (fabs(Bt39 - statrec.avg) > Bt39devavg)
	        flags_sst4[ip] |= SSTF_BTDEVAVG;
	    */
	}
        if (btboxstats(ip,ib40,NBANDSIR,btbox,btbox,&statrec) > 0) {
	    if ( (statrec.max - statrec.min) > Bt40unif1)
	        flags_sst4[ip] |= SSTF_BTNONUNIF;
	    if ( (statrec.max - statrec.min) > Bt40unif2)
	        flags_sst4[ip] |= SSTF_BTVNONUNIF;
	    /*
            if (fabs(Bt40 - statrec.avg) > Bt40devavg)
	        flags_sst4[ip] |= SSTF_BTDEVAVG;
	    */
	}

    }
    return;
}


/* ----------------------------------------------------------------------------------- */
/* set_quality_sst() - set quality levels for long wave sea surface temperature        */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
void set_qual_sst(l2str *l2rec)
{
    int32_t ip, ipb;
    float rhot; 

    for (ip=0; ip<l2rec->npix; ip++) {

        if (l2rec->solz[ip] < SOLZNIGHT) {

	    /* daytime 11um SST quality level */

            if ( 
                (flags_sst[ip] & SSTF_ISMASKED   ) > 0 ||
                (flags_sst[ip] & SSTF_BTBAD      ) > 0 )

                qual_sst[ip] = 4;

	    else if ( 
                (flags_sst[ip] & SSTF_VHISENZ    ) > 0 ||
                (flags_sst[ip] & SSTF_BTRANGE    ) > 0 ||
                (flags_sst[ip] & SSTF_SSTRANGE   ) > 0 ||
                (flags_sst[ip] & SSTF_BTVNONUNIF ) > 0 ||
	        (flags_sst[ip] & SSTF_SSTREFVDIFF) > 0 )

     	        qual_sst[ip] = 3;

            else if ( 
                (flags_sst[ip] & SSTF_BTNONUNIF ) > 0 ||
                (flags_sst[ip] & SSTF_REDNONUNIF) > 0 )

  	        qual_sst[ip] = 2;

            else if ( 
		(flags_sst[ip] & SSTF_SSTREFDIFF) > 0 ||
		(flags_sst[ip] & SSTF_HISENZ    ) > 0 )

  	        qual_sst[ip] = 1;

            else 

	        qual_sst[ip] = 0;

            /* Reduce quality if red-band reflectance is high and sst is cold (collect 5) */
	    if (haveRed) {
                ipb = ip*nbvis + ibred;
                rhot = PI * l2rec->Lt[ipb] / l2rec->Fo[ibred];
                if (qual_sst[ip] < 3 && sst[ip]-l2rec->sstref[ip] < -1.0 && rhot > 0.05)
	            qual_sst[ip]++;
	    }

	} else {

            if ( 
                (flags_sst[ip] & SSTF_ISMASKED   ) > 0 ||
                (flags_sst[ip] & SSTF_BTBAD      ) > 0 )

                qual_sst[ip] = 4;

	    else if ( 
                (flags_sst[ip] & SSTF_BTRANGE   ) > 0 ||
                (flags_sst[ip] & SSTF_SSTRANGE  ) > 0 ||
	        (flags_sst[ip] & SSTF_BT4REFDIFF) > 0 ||
	        (flags_sst[ip] & SSTF_SSTREFVDIFF) > 0 )

     	        qual_sst[ip] = 3;

            else if (  
                (flags_sst[ip] & SSTF_BTVNONUNIF) > 0 || 
	        (flags_sst[ip] & SSTF_SST4VDIFF ) > 0 ||
                (flags_sst[ip] & SSTF_VHISENZ   ) > 0 )

  	        qual_sst[ip] = 2;

            else if ( 
		(flags_sst[ip] & SSTF_SSTREFDIFF) > 0 ||
                (flags_sst[ip] & SSTF_BTNONUNIF ) > 0 ||   
	        (flags_sst[ip] & SSTF_SST4DIFF  ) > 0 ||
		(flags_sst[ip] & SSTF_HISENZ    ) > 0 )

  	        qual_sst[ip] = 1;

            else 

	        qual_sst[ip] = 0;

            /* reduce quality if 4um BTs show non-uniformity (collect 5) */
            if (qual_sst[ip] < 3 && (flags_sst4[ip] & SSTF_BTNONUNIF) > 0)
	        qual_sst[ip]++;


	}

    }
    return;
}


/* ----------------------------------------------------------------------------------- */
/* set_qual_sst4() - set quality levels for short wave sea surface temperature         */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
void set_qual_sst4(l2str *l2rec)
{
    int32_t ip;

    for (ip=0; ip<l2rec->npix; ip++) {

        if (l2rec->solz[ip] < SOLZNIGHT) {

	    /* daytime 4um SST quality level */

            if ( 
                (flags_sst4[ip] & SSTF_ISMASKED   ) > 0 ||
                (flags_sst4[ip] & SSTF_BTBAD      ) > 0 )

                qual_sst4[ip] = 4;

	    else    	    

                qual_sst4[ip] = 3;


	} else {

	    /* nighttime 4um SST quality level */

            if ( 
                (flags_sst4[ip] & SSTF_ISMASKED   ) > 0 ||
                (flags_sst4[ip] & SSTF_BTBAD      ) > 0 )

                qual_sst4[ip] = 4;

	    else if ( 
                (flags_sst4[ip] & SSTF_BTRANGE    ) > 0 ||
                (flags_sst4[ip] & SSTF_SSTRANGE   ) > 0 ||
                (flags_sst4[ip] & SSTF_BT4REFDIFF ) > 0 ||
	        (flags_sst4[ip] & SSTF_SSTREFVDIFF) > 0 )

     	        qual_sst4[ip] = 3;

            else if ( 
		(flags_sst4[ip] & SSTF_BTVNONUNIF) > 0 ||  
	        (flags_sst4[ip] & SSTF_SST4VDIFF ) > 0 ||
	        (flags_sst4[ip] & SSTF_VHISENZ   ) > 0 )

  	        qual_sst4[ip] = 2;

            else if ( 
		(flags_sst[ip]  & SSTF_SSTREFDIFF) > 0 ||
                (flags_sst4[ip] & SSTF_BTNONUNIF ) > 0 || 
	        (flags_sst4[ip] & SSTF_SST4DIFF  ) > 0 ||
		(flags_sst4[ip] & SSTF_HISENZ    ) > 0 )

  	        qual_sst4[ip] = 1;

            else 

	        qual_sst4[ip] = 0;
	}

    }
    return;
}


/* ----------------------------------------------------------------------------------- */
/* nlsst() - nonlinear sst, long wave sea surface temperature                          */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
float nlsst(float Bt11, float Bt12, float senz, float sstref, float lat, 
            float bounds[6][2], float a[6][4])
{
    static float pi = PI;
    static float dBtlo = 0.5;
    static float dBthi = 0.9;

    float dBt = Bt11 - Bt12; /* aka "thing" */
    float mu  = cos(senz*pi/180.0);
    float sstlo, ssthi, sst;
    int   ic;

    if ((evalmask & SSTMODS) != 0) {

	/* choose coeffs by latband */

	/* bounds go from -90 to 90 so we want the first set */
	/* that has max lat greater than pixel lat           */
	ic = 6;
	while (lat <= bounds[ic-1][1] + latwin) {
	    ic--;
	}
	dBtlo = bounds[ic][1]-latwin;
	dBthi = bounds[ic][1]+latwin;
	if (lat < bounds[ic][1]-latwin || ic == 5) {
	    sst = a[ic][0] + a[ic][1]*Bt11 + a[ic][2]*dBt*sstref + a[ic][3]*dBt*(1.0/mu-1.0);
	} else {
	    sstlo = a[ic][0] + a[ic][1]*Bt11 + a[ic][2]*dBt*sstref + a[ic][3]*dBt*(1.0/mu-1.0);
	    ssthi = a[ic+1][0] + a[ic+1][1]*Bt11 + a[ic+1][2]*dBt*sstref + a[ic+1][3]*dBt*(1.0/mu-1.0);
	    sst = sstlo + (lat-dBtlo)/(dBthi-dBtlo)*(ssthi-sstlo);
	}

    } else {
	if (dBt <= dBtlo)
	    sst = a[0][0] + a[0][1]*Bt11 + a[0][2]*dBt*sstref + a[0][3]*dBt*(1.0/mu-1.0);
	else if (dBt >= dBthi)
	    sst = a[1][0] + a[1][1]*Bt11 + a[1][2]*dBt*sstref + a[1][3]*dBt*(1.0/mu-1.0);
	else {
	    sstlo = a[0][0] + a[0][1]*Bt11 + a[0][2]*dBt*sstref + a[0][3]*dBt*(1.0/mu-1.0);
	    ssthi = a[1][0] + a[1][1]*Bt11 + a[1][2]*dBt*sstref + a[1][3]*dBt*(1.0/mu-1.0);
	    sst = sstlo + (dBt-dBtlo)/(dBthi-dBtlo)*(ssthi-sstlo);
	}
    }
 
    return(sst);
}

/* ----------------------------------------------------------------------------------- */
/* nlsst4() - nonlinear sst4, 4um sea surface temperature (3.9 & 4.0 um)               */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
float nlsst4(float Bt39, float Bt40, float senz, float lat, 
      float bounds[6][2], float a[6][5])
{
    int   ic;
    float dBt = Bt39 - Bt40;
    float mu  = cos(senz/RADEG);
    float sst4;
    float sst4lo;
    float sst4hi;
    float dBtlo;
    float dBthi;

    if ((evalmask & SSTMODS) != 0) {

	/* choose coeffs by latband */

	/* bounds go from -90 to 90 so we want the first set */
	/* that has max lat greater than pixel lat           */
	ic = 6;
	while (lat <= bounds[ic-1][1] + latwin) {
	    ic--;
	}
	dBtlo = bounds[ic][1]-latwin;
	dBthi = bounds[ic][1]+latwin;
	if (lat < bounds[ic][1]-latwin || ic == 5) {
	    sst4 = a[ic][0] + a[ic][1]*Bt39 + a[ic][2]*dBt + a[ic][3]*(1.0/mu-1.0) + a[ic][4]*dBt*(1.0/mu-1.0);
	} else {
	    sst4lo = a[ic][0] + a[ic][1]*Bt39 + a[ic][2]*dBt + a[ic][3]*(1.0/mu-1.0) + a[ic][4]*dBt*(1.0/mu-1.0);
	    sst4hi = a[ic+1][0] + a[ic+1][1]*Bt39 + a[ic+1][2]*dBt + a[ic+1][3]*(1.0/mu-1.0) + a[ic+1][4]*dBt*(1.0/mu-1.0);
	    sst4 = sst4lo + (lat-dBtlo)/(dBthi-dBtlo)*(sst4hi-sst4lo);
	}

    } else {
	sst4 = a[0][0] + a[0][1]*Bt39 + a[0][2]*dBt + a[0][3]*(1.0/mu-1.0) + a[0][4]*dBt*(1.0/mu-1.0);
    }
 
    return(sst4);
}


/* ----------------------------------------------------------------------------------- */
/* comp_sst4() - compute sea surface temperature                           .           */
/*                                                                                     */
/* B. Franz, SAIC, August 2003.                                                        */
/* ----------------------------------------------------------------------------------- */
void comp_sst4(l2str *l2rec)
{
    static int firstCall = 1;

    static float coef[6][5];
    static float bounds[6][2];

    int32_t  ip,ipb;
    float Bt39;
    float Bt40;
    float senz;
    float lat;
    float dBt;

    if (firstCall) {
        char *filedir;
        char filename[FILENAME_MAX];
        char sensor[20];
        if (l2rec->input->proc_sst != 1) {
            HY1C_out("-E- %s: SST processing must be enabled (proc_sst=1) to make sst4.\n", __FILE__);
            exit(1);
        }            



	/*
        if ((filedir = getenv("OCDATAROOT")) == NULL) {
            HY1C_out("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
            exit(1);
        }
        strcpy(sensor,&sensorName[l2rec->sensorID][0]);
        lowercase(sensor);
        strcpy(filename, filedir);
        if ((evalmask & SSTMODS) != 0) strcat(filename,"/eval");
        strcat(filename, "/");
        strcat(filename,sensor);
        strcat(filename,"/cal/");
        strcat(filename,"sst4_");
        strcat(filename,sensor);
        strcat(filename,".dat");
	*/

        if (l2rec->input->sst4coeffile != NULL) {
            read_sst4_coeff(l2rec->sensorID,l2rec->input->sst4coeffile,*(l2rec->year),*(l2rec->day),bounds,coef);
	} else {
	    HY1C_out("SST4 algorithm coefficient file not specified.");
            exit(1);
	}
        firstCall = 0;
    }


    for (ip=0; ip<l2rec->npix; ip++) {

        sst4[ip] = sstbad;
        flags_sst4[ip] = 0;

        /* skip if pixel already masked */
        if (sstmasked(l2rec,ip))
	    continue;

//        ipb = ip*NBANDSIR;
        ipb = ip*2;
        senz = l2rec->senz[ip];
        lat = l2rec->lat[ip];
        Bt39 = l2rec->Bt[ipb+ib39];
        Bt40 = l2rec->Bt[ipb+ib40];
        dBt  = Bt39-Bt40;

        /* skip pixel if BT cound not be computed */
        if (Bt39 < BT_LO+0.1 || Bt39 > BT_HI-0.1 || 
            Bt40 < BT_LO+0.1 || Bt40 > BT_HI-0.1) {
            continue;
        }

        /* compute SST4 */
        sst4[ip] = nlsst4(Bt39,Bt40,senz,lat,bounds,coef);
    }

    return; 
}


/* ----------------------------------------------------------------------------------- */
/* comp_sst() - compute sea surface temperature                            .           */
/*                                                                                     */
/* B. Franz, SAIC, August 2003.                                                        */
/* ----------------------------------------------------------------------------------- */
void comp_sst(l2str *l2rec)
{
    static int firstCall = 1;

    static float rvsint     = 0.01498041;   /* */
    static float rvsslope   = -0.001756443; /* */

    static float coef[6][4];
    static float bounds[6][2];
    static float dsst[2] = {0.0,0.0};

    float sstref;
    int32_t  ip, ipb;
    float Bt11;
    float Bt12;
    float senz;
    float solz;
    float lat;
    float dBt;
    float satzdir;
    float rvscor;
    int   mside = l2rec->mside;

    if (firstCall) {
        char *filedir;
        char filename[FILENAME_MAX];
        char sensor[20];
        if (l2rec->input->proc_sst != 1) {
            HY1C_out("-E- %s: SST processing must be enabled (proc_sst=1) to make sst.\n", __FILE__);
            exit(1);
        }            
        if (l2rec->input->sstcoeffile != NULL) {
            read_sst_coeff(l2rec->sensorID,l2rec->input->sstcoeffile,*(l2rec->year),*(l2rec->day),bounds,coef);
	} else {
	    HY1C_out("SST algorithm coefficient file not specified.");
            exit(1);
	}
        if (l2rec->input->sstmirrfile != NULL) {
            read_sst_mirror(l2rec->sensorID,l2rec->input->sstmirrfile,*(l2rec->year),*(l2rec->day),&dsst[0]);
	} 

        firstCall = 0;
    }
    

    for (ip=0; ip<l2rec->npix; ip++) {

        sst[ip] = sstbad;
        flags_sst[ip] = 0;

        /* skip if pixel already masked */
        if (sstmasked(l2rec,ip))
	    continue;

        ipb = ip*NBANDSIR;
        senz = l2rec->senz[ip];
        solz = l2rec->solz[ip];
        lat  = l2rec->lat [ip];
        Bt11 = l2rec->Bt[ipb+ib11];
        Bt12 = l2rec->Bt[ipb+ib12];
        dBt  = Bt11-Bt12;

        /* skip pixel if BT cound not be computed */
        if (Bt11 < BT_LO+0.1 || Bt11 > BT_HI-0.1 || 
            Bt12 < BT_LO+0.1 || Bt12 > BT_HI-0.1) {
            continue;
        }

        /* use 4um SST for night reference */
        if (haveSST4 && solz >= SOLZNIGHT && sst4[ip] != sstbad)
	    sstref = sst4[ip];
	else
            sstref = l2rec->sstref[ip];

	/* apply response vs scan angle correction (RVS) to AQUA sst */
	if ((evalmask && SSTMODS) != 0 && l2rec->sensorID == MODISA) {
   	    satzdir = (l2rec->pixnum[ip] < 677) ? -1.0 : 1.0;
	    rvscor = rvsint + (rvsslope * satzdir * senz);
	} else
	    rvscor = 0.0;
        /* compute SST & correct for mirror-side cooling differential */
        sst[ip] = nlsst(Bt11,Bt12,senz,sstref,lat,bounds,coef) + dsst[mside] + rvscor; 
    }

    return; 
}


/* ----------------------------------------------------------------------------------- */
/* run_sst() - run the SST algorithm over full scan and store in static arrays         */
/*                                                                                     */
/* Note: the 11um SST is using 4um SST as input (when we have it), so order matters.   */
/* Similarly, flags can't be computed until both SST products are computed.            */
/*                                                                                     */
/* B. Franz, SAIC, May 2005.                                                           */
/* ----------------------------------------------------------------------------------- */
void run_sst(l2str *l2rec)
{

    static int firstCall = 1;

    if (firstCall) {
        init_sst(l2rec);
        firstCall = 0;
    }  

    // see the line 1893
    // SSTFAIL/SSTWARN comment
    // 20180822
    return;

    /* compute SST */
    if (haveSST4)
        comp_sst4(l2rec);
    if (haveSST)
        comp_sst(l2rec);

    /* compute flags */
    if (haveSST4)
        set_flags_sst4(l2rec);
    if (haveSST)
        set_flags_sst(l2rec);

    /* set quality levels */
    if (haveSST4)
       set_qual_sst4(l2rec);
    if (haveSST)
        set_qual_sst(l2rec);

    /* set SSES */
    if (haveSSES) {
       if (haveSST4)
         set_sses_sst4(l2rec);
       if (haveSST)
         set_sses_sst(l2rec);
    }

    recnumSST = l2rec->iscan;
}




/* =================================================================================== */
/*              interface functions for l2_hdf_generic()                               */
/* =================================================================================== */

/* ----------------------------------------------------------------------------------- */
/* get_bias_sst() - SSES bias for long wave sea surface temperature                    */
/*                                                                                     */
/* B. Franz, SAIC, January 2006.                                                       */
/* ----------------------------------------------------------------------------------- */
float *get_bias_sst(l2str *l2rec)
{
    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    return(bias_sst);
}

/* ----------------------------------------------------------------------------------- */
/* get_bias_sst4() - SSES bias for short wave sea surface temperature                  */
/*                                                                                     */
/* B. Franz, SAIC, January 2006.                                                       */
/* ----------------------------------------------------------------------------------- */
float *get_bias_sst4(l2str *l2rec)
{
    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    return(bias_sst4);
}

/* ----------------------------------------------------------------------------------- */
/* get_stdv_sst() - SSES stdv for long wave sea surface temperature                    */
/*                                                                                     */
/* B. Franz, SAIC, January 2006.                                                       */
/* ----------------------------------------------------------------------------------- */
float *get_stdv_sst(l2str *l2rec)
{
    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    return(stdv_sst);
}

/* ----------------------------------------------------------------------------------- */
/* get_stdv_sst4() - SSES stdv for short wave sea surface temperature                  */
/*                                                                                     */
/* B. Franz, SAIC, January 2006.                                                       */
/* ----------------------------------------------------------------------------------- */
float *get_stdv_sst4(l2str *l2rec)
{
    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    return(stdv_sst4);
}

/* ----------------------------------------------------------------------------------- */
/* get_qual_sst() - quality flags for long wave sea surface temperature                */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
int8 *get_qual_sst(l2str *l2rec)
{
    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    return(qual_sst);
}

/* ----------------------------------------------------------------------------------- */
/* get_qual_sst4() - quality levels for short wave sea surface temperature             */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
int8 *get_qual_sst4(l2str *l2rec)
{
    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    if (!haveSST4) {
        HY1C_out("short-wave SST is not available for this sensor.\n");
        exit(1);
    }

    return(qual_sst4);
}

/* ----------------------------------------------------------------------------------- */
/* get_flags_sst() - quality flags for long wave sea surface temperature               */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
int16 *get_flags_sst(l2str *l2rec)
{
    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    return(flags_sst);
}

/* ----------------------------------------------------------------------------------- */
/* get_flags_sst4() - quality flags for short wave sea surface temperature             */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
int16 *get_flags_sst4(l2str *l2rec)
{
    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    if (!haveSST4) {
        HY1C_out("short-wave SST is not available for this sensor.\n");
        exit(1);
    }

    return(flags_sst4);
}

/* ----------------------------------------------------------------------------------- */
/* get_sst() - return longwave sea surface temperature for current scan                */
/*                                                                                     */
/* B. Franz, SAIC, August 2005.                                                        */
/* ----------------------------------------------------------------------------------- */
float *get_sst(l2str *l2rec)
{
    int32_t  ip;

    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    if (!haveSST) {
        HY1C_out("long-wave SST is not available for this sensor.\n");
        exit(1);
    }

    for (ip=0; ip<l2rec->npix; ip++) {
        switch (qual_sst[ip]) {
	  case 4:
            // 20180822
            // l2rec->flags[ip] |= SSTFAIL;
            break;
	  case 3:
	  case 2:
	  case 1:
            //20180822
            // l2rec->flags[ip] |= SSTWARN;
            break;
	}
    }

    return(sst);
}

/* ----------------------------------------------------------------------------------- */
/* get_sst4() - return shortwave sea surface temperature for current scan              */
/*                                                                                     */
/* B. Franz, SAIC, May 2005.                                                           */
/* ----------------------------------------------------------------------------------- */
float *get_sst4(l2str *l2rec)
{
    int32_t  ip;

    if (!sst_ran(l2rec->iscan))
        run_sst(l2rec);

    if (!haveSST4) {
        HY1C_out("short-wave SST is not available for this sensor.\n");
        exit(1);
    }

    for (ip=0; ip<l2rec->npix; ip++) {
        switch (qual_sst4[ip]) {
	  case 4:
            l2rec->flags[ip] |= SSTFAIL;
            break;
	  case 3:
	  case 2:
	  case 1:
            l2rec->flags[ip] |= SSTWARN;
            break;
	}
    }

    return(sst4);
}

